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We calculate the antideuteron yield in dark matter annihilations on an event-by-event basis using 
the HERWIG+- 1- Monte Carlo generator. We present the resulting antideuteron fluxes for quark and 
gauge boson final states. As deuteron production in the coalescence model depends on momentum 
differences between nucleons that are small compared to Aqcd, it is potentially very sensitive to 
the hadronization model employed. We therefore compare our antideuteron yields to earlier results 
based on PYTHIA, thereby estimating their uncertainties. We also briefly discuss the importance 
of n > 2 final states for annihilations of heavy DM particles. 
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I. INTRODUCTION 

Despite various cosmological and astrophysical indica- 
tions for the presence of nonbaryonic dark matter (DM) , 
its particle nature has yet to be proven. Restricting the 
space of candidates to thermal relics, the thermally aver- 
aged annihilation cross section (av) at freeze-out is fixed 
by the DM abundance, while the DM mass Mdm is only 
weakly constrained: Unitarity of the S'-matrix constrains 
the mass of any thermal relic as m < 50 TeV [T] , while the 
requirement that the DM is cold translates for thermal 
relics into a lower mass limit of the order 10 keV. 

One strategy towards DM detection is to search for 
its self-annihilation (or decay) products. The annihila- 
tion of symmetric DM leads to equal injection rates of 
matter and antimatter particles into the Galaxy, while 
the cosmic ray flux from astrophysical sources is matter- 
dominated. A possible way to detect DM is therefore to 
carefully estimate the expected antimatter fluxes from 
astrophysical sources, and then to search for any excess. 
The authors of Ref. [2] suggested antideuterons as a sig- 
nature in addition to the usually discussed antiproton 
and positron signal from DM annihilations. In particu- 
lar, they argued that the DM antideuteron spectra are 
much flatter at low energies than those from cosmic-ray- 
gas interactions. 

The correct application of the coalescence model in 
deuteron production requires its implementation on an 
event-by-event basis 01!]. Since deuteron production de- 
pends on momentum differences po that are small or com- 
parable to Aqcd ~ 200 MeV, one may expect a rather 
strong dependence of the results on the hadronization 
models employed in the Monte Carlo simulation. More- 
over, the coalescence model probes not only the energy 
spectrum of nucleons, but also their two-particle correla- 
tions in momentum space. Since the physics underlying 
different hadronization models — e.g. cluster hadroniza- 
tion versus string fragmentation — varies strongly, the 
choice of hadronization model could thus have a profound 
effect on the generated antideuteron spectra. 

The main aim of this work is therefore to derive re- 
sults for the antideuteron yields calculated using HER- 



WIG++ [5] version 2.4.2 (based on a cluster hadroniza- 
tion model) and to compare them to earlier results using 
PYTHIA [6] (based on string fragmentation). We also 
briefly discuss the importance of n > 2 final states for 
annihilations of heavy DM particles. 

II. DEUTERON PRODUCTION IN THE 
COALESCENCE MODEL 

The fusion of (anti)protons and (anti)neutrons into 
(anti)deuterons is usually described with the so-called co- 
alescence model. This model is based on the assumption 
that any nucleons with a momentum difference Ap < po 
for some given p$ will merge and form a nucleus. In 
its initial form, the model was applied to deuteron pro- 
duction in nucleus-nucleus interactions. In the simplest 
approximation, the final-state hadrons in nucleus-nucleus 
scattering are formed in a "fireball" and are emitted close 
to isotropically in the center-of-mass (CoM) frame. As- 
suming additionally that correlations between nucleons 
are negligible, the antideuteron energy spectra in the lab 
frame can be derived from the energy spectra of nucleons 
as [71 E| 

dN s = pljmd 1 dN n dN p 

dT 6 m n m p + 2mdT(l dT dT ' 

Here, po is the maximal momentum difference allowed 
between a pn pair in order to form an antideuteron ac- 
cording the coalescence model, m^, m p , and m n denote 
the deuteron, proton and neutron masses respectively, 
and the nucleon spectra on the RHS are to be evalu- 
ated at the value Tp = T fl — Tj/2 of the kinetic energy 
Ti = Ei- m.j. 

The constant prefactor in Eq. varies depending on 
the assumptions made in its derivation [J] . Such a factor 
can, however, be absorbed by re-defining pq. For later 
reference, we note that our definition agrees with the one 
of Ref. 3], while the one of Refs. [SUTU] differs by a factor- 
eight. 

The coalescence model can easily be implemented di- 
rectly in a Monte Carlo simulation by comparing the mo- 
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menta of the final state nucleons in their respective CoM 
frames for each individual annihilation event. We will 
from now on refer to results from calculations where co- 
alescence was applied on an event-by-event basis within 
the Monte Carlo as "Monte Carlo" results, while results 
from calculations where coalescence was applied to the 
average antiproton and antineutron energy spectra using 
Eq. ([IJ will be referred to as "isotropic" results. Corre- 
spondingly, we refer to the two approaches as the Monte 
Carlo approach and the isotropic approach. 



A. Determining po 

The momentum threshold po for the formation of 
deuterons is found in both approaches by running sim- 
ulations and adjusting po such that the computational 
result matches the experimental one. Antideuteron pro- 
duction in e + e~ -collisions was studied by the ALEPH 
collaboration [TT] using LEP-I data. At the Z-resonance, 
each hadronic event was found to give rise to (5.9 ± 
1.8 ± 0.5) x 10~~ 6 antideuterons in the momentum range 
0.62 GeV < % < 1.03 GeV in the angular acceptance 
range | cosi?| < 0.95 of the detector. 

Our simulations using HERWIG++ reproduce the ex- 
perimentally measured antideuteron yield choosing p = 
110 MeV in the "Monte Carlo approach" and po = 
126 MeV in the "isotropic approach", respectively. Note 
that the po values differ in general for the two methods, 
and should therefore be self-consistently calibrated. 

There is a significant range in the value of po used in 
the literature. Reference [9] suggests a window between 
66 MeV and 105 MeV for the isotropic approach (note 
that their definition of p requires a rescaling of p by 
2 to compare with our values), while Kadastik et al. [3] 
use a value of po = 160 MeV for both the isotropic and 
Monte Carlo approach. 

The relatively large difference of 45% between the po 
values found by us using HERWIG++ and Kadastik et 
al. using PYTHIA indicates that the differing hadroniza- 
tion models do indeed lead to physically different predic- 
tions. Since the momentum range used for the calibra- 
tion is very small, our normalization procedure does not 
even guarantee that the total number of antideuterons 
predicted by HERWIG++ and PYTHIA for our refer- 
ence process e + e~ — > Z° — > hadrons agrees. In order 
to illustrate this point, we repeated the calibration us- 
ing PYTHIA (version 8.160). In Fig. [IJ we compare the 
energy spectrum TdN/dT of antideuterons calculated by 
us with the two QCD simulations. By construction, the 
two spectra cross in the calibration range, but PYTHIA 
predicts a harder energy spectrum and a total yield of 
antideuterons larger by a factor two than HERWIG. Go- 
ing to higher CoM energies or DM masses, one expects 
moreover that the spectral shape and the total number of 
antideuterons evolves differently in the two hadronization 
schemes. 

Note that the single nucleon spectra calculated with 
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FIG. 1: Comparison of the energy spectrum TdN/dT of an- 
tideuterons for the calibration reaction, with and without the 
angular cut cos i9 < 0.95. The energy range used for the 
calibration is shown as the green area. 

HERWIG ++ and PYTHIA agree quite well, implying 
that the large differences found in the antideuteron yield 
are caused by differences in the two-particle correlations 
of the ph pairs. 



B. Fragmentation spectra 

Our results were generated using 10 s events for each 
annihilation channel and dark matter mass, and are gen- 
erally presented in terms of the scaled kinetic energy 
x = T/MdMi where Mdm is the dark matter mass. 
For practical reasons, we used the lightest neutralino 
in the MSSM as DM particle, thus allowing us to use 
MadGraph [TJ] to generate events to be fed into HER- 
WIG++ for showering and hadronization. Note that 
while the branching ratios into different final states are 
strongly model-dependent, the energy spectrum of a spe- 
cific final-state is — except in cases where spin correlations 
are important — mainly determined by the CoM energy 
y/8 = 2M DM . 



1. The antideuteron fragmentation spectra 

The antideuteron fluxes from the isotropic and Monte 
Carlo approaches obtained by us using HERWIG++ arc 
plotted for DM masses of 300 GeV and 1 TeV in Fig. [2] In 
the bb case, the low- a: part of the energy spectra dN/dx 
are roughly the same for the two methods. The spectra 
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in the Monte Carlo approach extend, however, towards 
much larger x. For the W + W~ case, there is a significant 
difference in the antideuteron spectrum for all x between 
the two approaches. For 1 TeV, the difference is of order 
100 at low x, while it is around a factor ~ 30 for 300 GeV. 
The shapes of the spectra also differ somewhat between 
the two approaches, with the maxima of xdN/dx shifted 
towards higher x in the Monte Carlo approach. 

It is clear that when using the correct Monte Carlo 
approach, the difference in magnitude between the bb 
and W + W~ antideuteron spectra becomes much smaller 
than in the isotropic approach. Another interesting fea- 
ture is that this difference appears to depend on the DM 
mass. 

In order to investigate the dependence of the an- 
tideuteron spectra on the DM mass, we plot the to- 
tal number of antideuterons produced against the DM 
mass for the two approaches in Fig. [3j Comparing this 
graph to the corresponding one of Ref. [3] obtained us- 
ing PYTHIA, we find that the general behavior of the 
antideuteron yields as function of channel and approach 
agrees well. Comparing individual energy spectra, we ob- 
serve a sharper drop at high energies in the spectra pro- 
duced by HERWIG++. We also note some differences 
in the overall magnitudes, and we therefore proceed by 
comparing our HERWIG++ results to those of Kadastik 
et al. [5] in more detail. For their results we use the data 
files (without EW corrections) from Ref. [H] , which were 
generated using PYTHIA version 8.135. 

In Fig. |4j we show the ratio 

dN/dx\ HERWIG++ 

H = 



bh source spectrum 



dN/dx\ 



PYTHIA 



of the antideuteron energy spectra calculated with HER- 
WIG++ and PYTHIA, respectively. In the peak region, 
this ratio is close to 0.5 for the W + W~ channel, while 
it is close to 0.3 for the bb channel. In the forward re- 
gion x — > 1, PYTHIA predicts a significantly higher an- 
tideuteron yield than HERWIG++, leading a drop of 
the ratio R. As expected from the calibration reaction, 
cf. Fig. [T] the ratio R lies below one in the peak region, 
followed by a crossover to R > 1 at lower x. The ex- 
pected crossover is not seen in the 300 GeV quark case, 
possibly due to insufficient statistics at low x. For the 
gauge bosons, the x value for which the ratio crosses one 
appears to be shifted slightly toward higher values with 
increasing DM mass, but the uncertainty with 10 s events 
is too large to determine the crossover points precisely. 
For the quark case, the crossing region is shifted towards 
lower x. In order to see if there is a DM mass dependence 
of the crossing point in the quark channel, we also inves- 
tigate the ratio at A/dm = 50 GeV, as seen in Fig. [5| In 
addition to the 6-quarks studied in the rest of the article, 
the figure also includes the corresponding ratio for anni- 
hilations into light quarks (u, d, s). We see that not only 
does the crossover for the 6-quark case appear at higher x 
for lower A/dmj but the ratios of the antideuteron energy 
spectra in general also differ somewhat between light and 
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FIG. 2: Antideuteron source spectra for dark matter annihi- 
lations into bb (top) and W + W~ (bottom). The solid lines 
show the spectra for per-event coalescence within the Monte 
Carlo, while the dotted lines show the spectra for coalescence 
of the average antiproton and antineutron spectra. Red lines 
show the result for a dark matter mass of 1 TeV, blue lines 
for 300 GeV. 



heavy quarks. 



2. Interpretation 



The perhaps most striking feature of Fig. [3] is how 
strongly the impact of going from the isotropic to the 
Monte Carlo approach varies between the bb and W + W~ 
channel. In the latter case, changing the DM mass by 
a factor 10 leads to a factor 100 difference in the an- 
tideuteron yield predicted by the two approaches, while 
in the former case the prediction differs only by a factor 
three. This suggests that the basic assumption under- 
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FIG. 3: Average total antideuteron yield per DM annihila- 
tion event. The solid lines show the results using per-event 
coalescence, while the dashed lines show the results using the 
isotropic approximation. 
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FIG. 4: Ratio R of the antideuteron spectra predicted by 
HERWIG++ and PYTHIA as function of x for neutralino 
annihilations with mass 300 GeV and 1 TeV, respectively. 



lying the "isotropic approach," namely an uncorrelated 
isotropic emission of the final-state particles, is not too 
strongly violated for qq final-states. 

To illustrate this effect, we show in Fig. [6] the distri- 
bution of the angles between the momenta of antiproton 
and antineutron pairs. These distributions are strongly 
peaked for angles near and tt for the W + W~ channel, 
while the peaks are much less pronounced for bb. More 
importantly, in the quark channel the shape of the dis- 
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FIG. 5: Ratio R of the antideuteron spectra predicted by 
HERWIG++ and PYTHIA as function of x for X X -> bb and 
XX — ^ with m x = 50 GeV and q = u,d, s. 



tribution changes only mildly when increasing the DM 
mass, while the peaks in the W + W~~ channel increase 
strongly. 

The latter effect is a simple reflection of relativistic 
beaming: The W's fed into HERWIG++ from Mad- 
Graph are decaying on-shell, and their decay product are 
emitted in a cone of opening angle $ ~ mw/MnM- The 
relativistic beaming in itself does not affect coalescence 
in the Monte Carlo approach, as the nucleon momenta 
are to be evaluated in the CoM frame of the respective 
pn pairs. The fact that the bosons are decaying on-shell, 
however, implies that the nucleon multiplicity — and thus 
the antideuteron multiplicity — remains constant with in- 
creasing Mdm, as shown in Fig. [7j Using the (wrong) 
Eq. ([I]) derived in the isotropic approach, 



dx 



1 dNft dN* 



oc 



Km 



dx dx 



(2) 



one would instead expect an f/M^ M suppression in 
the antideuteron yield, as is reflected by the isotropic 
W+W~ result in Fig. || 

In contrast to the gauge bosons, the &'s are treated 
as virtual particles and start QCD cascades. Although 
these cascades are angular-ordered and lead to confined 
jets, the overall shape of the event is not too far from 
spherically symmetric. Moreover, the two initial back- 
to-back jets are color connected, and partons from the 
two jets have to combine in order to produce colorless 
final states. 

The fact that the antideuteron yield does not show the 
I/Mp M suppression in the bb channel using the isotropic 
approach can also be explained. In this case, the growth 
of the nucleon multiplicity shown in Fig. [7] with Mtjm 
overcompensates the l/Afp M factor, leading to a slight 



5 



increase of the antideuteron yield even in the isotropic 
approach. 

Finally, we discuss the ratios of the HERWIG++ and 
PYTHIA fragmentation spectra. As already mentioned, 
the x for which R crosses one in the W + W~ channel ap- 
pears to shift slightly towards higher x with increasing 
DM mass: As the gauge bosons decay on-shell, their de- 
cay products are boosted according to the DM mass, and 
for Mdm 3> tnw we expect Xj to increase linearly with 
Mdm- With the T- value of the crossover point increas- 
ing linearly with Mdm, the corresponding x = T/A/dm 
should then converge towards a constant value with in- 
creasing DM masses. 

In the bb channel, however, we see from Figs. ||and |] 
that the ir-value of the crossing region appears to de- 
crease as ~ 1/MdMj such that the crossing occurs at a 
roughly constant value of the antideuteron kinetic energy 
Tj. Note, however, that the drop in the ratio for x — >• 1 
appears at a constant value of x rather than T. 

Using the differences between HERWIG++ and 
PYTHIA as an estimate for the uncertainty introduced 
by the hadronization models, we conclude that the high- 
energy part of the antideuteron spectra cannot be pre- 
dicted reliably for neither the gauge boson nor quark 
channels. For the bb channel, the spectrum has an 
uncertainty of a factor ~ 3 from x ~ 10 _1 down to 
T ~ 10~ 2 GeV, below which the uncertainty increases. 
For the gauge bosons, we find an uncertainty of ~ 2 
for 10 -3 < x < 0.5, with rapidly increasing uncertain- 
ties outside this range. For the W + W~ channel, the 
crossover point, and consequently the high uncertainty 
region at low energies, moves linearly towards higher T 
with increasing DM mass. Since T, rather than x is 
the relevant quantity for the observable intensity of an- 
tideuterons, this implies that observationally relevant en- 
ergies could lie in the region with large uncertainties for 
DM masses in the TeV range and above. 



C. Higher order processes 

Motivated by the PAMELA excess, previous works 
have discussed the annihilation and the resulting an- 
tideuteron fluxes of DM particles up to 30 TeV [3 [TO]. 
The unitarity limit of 50 TeV implies that the partial 
wave amplitude for the annihilation of such a heavy DM 
particles is close to one. Thus DM particles with masses 
> 10 TeV should behave as strongly interacting particles, 
and higher order processes should give an important con- 
tribution to their total annihilation cross section. 

As a test of this conjecture, we calculated the anni- 
hilation cross sections for XiXi ~~ ► W + W~, and for the 
higher order tree-level processes in which one and two 
extra Z-bosons are emitted. The calculations were per- 
formed for several different neutralino masses in the range 
200 GeV to 2 TeV. Using these cross sections, we calcu- 
lated the branching ratios for these processes, normalized 
to a sum of 1. The results from these calculations are 



plotted in Fig. [Sj 

Figure [8] shows as expected that the relative contribu- 
tion from higher order processes is negligible for low neu- 
tralino masses. The contributions do, however, increase 
rapidly with increasing masses. For 2 TeV, the branch- 
ing ratio for W + W~ Z is roughly 10% of that for the tree 
level process. While the process involving two additional 
Z-bosons is more strongly suppressed for low neutralino 
masses than the single Z-boson process, its relative con- 
tribution increases faster with increasing masses than for 

w+w-z. 

Extrapolating our results, we expect the contribution 
from higher order processes to become significant for a 
neutralino mass in the 10 TeV range [15]. The emission of 
additional Z-bosons are, of course, not the only possible 
higher order processes, and when performing calculations 
in the multi TeV range, the contributions from a number 
of processes should be considered. 

Note also that the results of this subsection are more 
mo del- dependent than the previous ones. For DM models 
other than the MSSM, the DM mass where final-states 
with n > 2 particles become important may be lower. 

III. ANTIDEUTERON INTENSITY AT EARTH 

We now move from the antideuteron spectra produced 
in a single annihilation to the calculation of the result- 
ing antideuteron intensity at Earth. As charged parti- 
cles scatter on fluctuations of the Galactic magnetic field 
with variation scales comparable to their Larmor radius, 
their propagation resembles a random walk and is well 
described by the diffusion approximation. We model this 
random walk using the so-called two-zone propagation 
model for the Galaxy. This is a cylindrical diffusion 
model consisting of a magnetic halo and a thin gaseous 
disk of radius R = 20 kpc and half-heights L (a free pa- 
rameter) and h = 100 pc, respectively. For the numerical 
values of the parameters in this model, we adopt the three 
sets presented in |16) to yield maximal, median and mini- 
mal antiproton fluxes from DM annihilations while being 
compatible with the observed B/C ratio. These sets are 
labeled 'max', 'med', and 'min', respectively, and listed in 
Table [TJ Assuming steady state conditions, the diffusion 
equation neglecting reacceleration is given by 

- D(T)V 2 f + |-(sign(z)/U c ) =Q- 2hS(z)T ann (T)f , 
oz 

(3) 

where f(x,T) = dNg/dT is the number density of an- 
tideuterons per unit kinetic energy, D(T) = D (31Z S the 
(spatial) diffusion coefficient, V c a convective wind per- 
pendicular to the Galactic disk, z the vertical coordinate, 
P — v/c the velocity and 1Z the rigidity of antideuterons. 

The term r ann (T) in Eq. ^ accounts for annihilations 
of antideuterons on the interstellar medium. The anni- 
hilation rate is given by 

r a nn = (™ff+4l?l H c)(<>>, (4) 
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FIG. 6: Distribution of the angles between the momenta of antiprotons and antineutrons; left for A/dm = 300 GeV, right for 
A/dm = 1 TeV. Blue and green lines show the distribution in "Monte Carlo" approach result for bb and W + W~ , respectively, 
while the red lines show the result of an isotropic distribution for comparison. 
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FIG. 7: Average total antiproton and antineutron yields per 
annihilation event as function of the dark matter mass. Green 
indicates antineutrons, while blue indicates antiprotons. The 
solid lines show the results for the bb case, while the dashed 
lines show the results for the W + W~ case. 



where the factor 4 a accounts for the different cross sec- 
tions of H and He interactions assuming simple geometri- 
cal scaling, and we use i*h ~ 1 cm -3 and 71h ~ 0.07tih as 
the number density of hydrogen and helium in the disk, 
respectively. The fit to the experimental data [13 US] we 
used for these reactions is shown in Fig. [9] 

Finally, the source term Q is fixed by the DM halo 



FIG. 8: Branching ratios for various annihilation channels 
as function of the dark matter (neutralino) mass. The red 
line shows the benchmark branching ratio for the tree level 
process XiXi ~^ W + W~ . The blue and green lines show 
the branching ratios for the higher order processes XiXi ~~ ► 
W + W~ZZ, respectively. 
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profile, which we choose as either a Navarro-Frenk- White 
(NFW) profile 17 (a = 1, /3 = 3, 7 = 1) or an isothermal 
profile (a = 2, /3 = 2, 7 = 0) in 
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max 


15 


0.46 
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med 
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0.0112 
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min 


1 


0.85 
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TABLE I: Propagation parameters for the max, med and min 
models. 



DM Halo Profile 


Po in GeV/cm 3 


a in kpc 


NFW 


0.26 


20 


Isothermal 


1.16 


5 


Einasto 


0.06 


20 



TABLE II: Density profile parameters. 



or as the Einasto profile, 



PEinastoM = PO exp 



2,,ry 
a \\a/ 



- 1 



a = 0.17. (6) 



Values of the free parameters po and a in these profiles 
are given for the Milky Way in Table |llj 

The solution to the diffusion equation Eq. ^ gives the 
intensity Jj at the position of the Earth as [Ml [T5] 



B 



f-?-V 



R(T) 



2 dT 



(7) 



where B is a possible boost factor, and the astrophysics 
is encoded in the propagation function R(T). Our results 
for R(T) are shown in Fig. [To] and agree well with those 
of Kadastik et. al. |3] (when exchanging the labels for 
the Einasto and Moore profiles in their figure). 

Combining the propagation function R(T) with the 
calculated energy spectra of antideuterons and account- 
ing for solar modulations by replacing T with T Q = 
T— |Ze|(/>Fisk with <^Fisk = 0.5 GV [5D], we are in position 
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FIG. 9: Cross section data for antideuterons on interstellar 
protons as a function of the antideuteron momentum. The 
points indicate experimental data, while the lines show the 
fits to the data which were used in our calculations. 
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FIG. 10: The function R(T), plotted for different dark matter 
profiles and propagation settings. The filled grey areas show 
the differences in R(T) between the density profiles for a given 
propagation model. The upper lines correspond to the 'max' 
model, the middle lines correspond to the 'med' model, and 
the lower lines correspond to the 'min' model. 



to determine the intensity of antideuterons at the posi- 
tion of the Earth. For the final results, we adopt for the 
thermally averaged cross section (av) = 3 x 10 _26 cm 3 /s, 
a boost factor B = 1 for each annihilation channel, and 
use the NFW density profile. 

The resulting intensity of antideuterons using the 
'med' propagation parameters is shown in Fig. |11| We 
see that there is a significant enhancement in the peaks of 
the spectra going from the isotropic to the correct Monte 
Carlo approach. This enhancement is most significant 
for the 1 TeV W + W~ case, where the peak in the Monte 
Carlo approach is two orders of magnitudes higher than 
in the isotropic approach. 

In Fig. 12 we compare the intensities predicted by 
HERWIG++ and PYTHIA simulations, giving an es- 
timate on the uncertainty in the hadronization. From 
the figure, we see that the uncertainty in the propaga- 



tion model leads to a spread in the final antideuteron 
spectrum of ~ 1.5 orders of magnitude, while the uncer- 
tainty in hadronization generally leads to an uncertainty 
of a factor 2-4, depending on the process in question. 
For most energies, the uncertainty in the propagation 
clearly dominates over the uncertainty from hadroniza- 
tion. For the quark channel, however, the uncertainty 
in hadronization becomes competitive for energies corre- 
sponding to x > 10 _1 . Because of the shift in kinetic 
energy due to solar modulation, the sharp rise in uncer- 
tainty at low energies in the W + W~ channel is not visible 
in Fig. 12 even with extended plot ranges. We expect, 
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however, that this uncertainty will become important for 
relevant ranges of the spectrum for DM masses in the 
TeV range. 

IV. CONCLUSIONS 

We calculated the antideuteron yield in dark mat- 
ter annihilations in the coalescence model on an event- 
by-event basis and presented the resulting antideuteron 
fluxes for quark and gauge boson final states. We 
showed that deuteron production is very sensitive to the 
hadronization model employed. Comparing our results 
using the HERWIG++ Monte Carlo simulation to ear- 
lier results using PYTHIA, we found that the predicted 
antideuteron yield varies by a factor ~ 3 for x < 0.1, 
T > 0.01 GeV the quark channel, and a factor ~ 2 in 
the range 10~ 3 < x < 0.5 in the gauge boson channel. 
Outside these ranges, the uncertainty increases rapidly 



in both channels. Using the differences between HER- 
WIG++ and PYTHIA as an estimate for the uncertainty 
introduced by the hadronization models, this results in an 
uncertainty of a factor of 2-4 on the observable intensity. 
While this is generally subdominant to the uncertainty 
from the propagation model at around 1.5 orders of mag- 
nitude, it is potentially dominant for the gauge boson 
channel at low energies for DM masses in the TeV range 
and above. We also showed the importance of n > 2 fi- 
nal states for annihilations of DM particles with masses 
> 10 TeV. 
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FIG. 11: Final antideuteron spectra near Earth after propagation and Solar modulation. Calculations are done for dark matter 
masses of 1 TeV and 100 GeV, using the NFW density profile and the 'med' propagation parameters. In both plots, we assumed 
a thermally averaged cross section of (av) = 3 x 10 _26 cm 3 /s. Continuous lines show the result for the Monte Carlo approach, 
while dashed lines show the result from the isotropic approach. 
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FIG. 12: Comparison of the antideuteron spectra from HERWIG++ and PYTHIA near Earth after propagation and Solar 
modulation. The plots show the results for dark matter masses of 1 TeV and 300 GeV, using the NFW density profile and 
varying propagation parameters. A thermally averaged cross section of (av) = 3 x 10~ 2e cm 3 /s. was assumed. 



